Available online at www.sciencedirect.com 


ScienceDirect 


JOURNAL OF 


www.elsevier.com /locate /jpowsour 


ELSEVIER 


Journal of Power Sources 173 (2007) 346-356 


AC impedance modelling study on porous electrodes of proton exchange 
membrane fuel cells using an agglomerate model 


Dietmar Gerteisen*, Alex Hakenjos, Jiirgen O. Schumacher 


Fraunhofer Institute for Solar Energy Systems, Heidenhofstr. 2, 79110 Freiburg, Germany 


Received 16 March 2006; received in revised form 15 March 2007; accepted 29 April 2007 
Available online 5 May 2007 


Abstract 


A one-dimensional model of the PEM fuel cell cathode is developed to analyse ac impedance spectra and polarisation curves. The porous gas 
diffusion electrode is assumed to consist of a network of dispersed catalyst (Pt/C) forming spherically shaped agglomerated zones that are filled 
with electrolyte. The coupled differential equation system describes: ternary gas diffusion in the backing (O2, N2, water vapour), Fickian diffusion 
and Tafel kinetics for the oxygen reduction reaction (ORR) inside the agglomerates, proton migration with ohmic losses and double-layer charging 
in the electrode. Measurements are made of a temperature-controlled fuel cell with a geometric area of 1.4cm x 1.4cm. Lateral homogeneity is 
ensured by using a high stoichiometry of Amin. The model predicts the behaviour of measured polarisation curves and impedance spectra. It is 
found that a better humidification of the electrode leads to a higher volumetric double-layer capacity. The catalyst layer resistance shows the same 
behaviour depending on the humidification as the membrane resistance. Model parameters, e.g. Tafel slope, ionic resistance and agglomerate radius 
are varied. A sensitivity analysis of the model parameters is conducted. 
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1. Introduction 


A proton exchange membrane (PEM) fuel cell is an elec- 
trochemical device in which the energy of a chemical reaction 
is converted directly into electricity. Hydrogen fuel is used in 
combination with air to produce electricity without combus- 
tion. The only by-products of the electrochemical reactions in 
the fuel cell are heat and water. Fuel cells are clean, quiet and 
energy-efficient. Possible application fields that are currently 
considered for PEM fuel cells include portable applications, sta- 
tionary power generation and the automotive area. PEM fuel 
cells use a solid polymer membrane as electrolyte. The electro- 
chemical reactions take place in porous gas diffusion electrodes. 

Several processes are involved during fuel cell operation: the 
electrochemical reactions in the electrodes, transport of elec- 
trons and protons, mass-transport of the gas species and water, 
and heat transport. Electrochemical impedance spectroscopy is 
a powerful characterisation method to understand these com- 
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plex processes in the fuel cell and electrodes, respectively. This 
measurement technique is based on the fact that the processes 
mentioned above take place on different characteristic time 
scales. Applying a small alternating voltage signal on a dc load 
level, the current response characterised by the amplitude and 
phase shift depends on the relaxation time of the coupled pro- 
cesses and thus on the perturbation frequency. The complex 
impedance is defined by the quotient of the perturbation voltage 
input and the ac current response. The impedance at different 
frequencies is called an impedance spectrum. The effects of the 
reaction kinetics, mass-transport processes and the limited ionic 
conductivity through the catalyst layer (CL) can be separated 
by analysing measured impedance spectra with a mathematical 
model. 

Different approaches can be found in the literature to describe 
the dc and ac characteristics of a PEM fuel cell cathode. Math- 
ematical models based on differential equations that describe 
the physics and the electrochemistry differ mostly in their geo- 
metrical properties. Springer et al. modeled the cathode as a 
homogeneous electrode and calculated polarisation curves [1]. 
Moreover, they presented a homogeneous cathode model for 
the calculation of ac impedance spectra [2]. Giner and Hunter 
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Nomenclature 

b Tafel slope (V) 

c total gas concentration in the backing (mol m~*) 

c* concentration of the dissolved oxygen at the 
agglomerate surface (mol m~?) 

c/c* normalised local oxygen concentration 

CpL,v volumetric double-layer capacity (F m~?) 

D diffusion coefficient of dissolved oxygen in an 
agglomerate (m? s~!) 

f Laplace-transformed function f 

F Faraday constant (9.6485 x 10* C mol!) 

ig exchange current density of the ORR at partial 
pressure of 1 atm (A m~? atm7!) 

i ig I”, the product of the exchange current density 
of the ORR at partial pressure of | atm with the 
ratio of the active surface area inside an agglomer- 
ate per unit agglomerate volume (A m~? atm!) 

Jcell fuel cell current density (A m~?) 

Jo, oxygen flux inside an agglomerate (mol m~? s~!) 

Jp protonic current density inside an agglomerate 
(Am?) 

J protonic current density in the catalyst layer 

L thickness of the catalyst layer (m) 

n number of transferred electrons 

POs oxygen partial pressure (atm) 

Pa operating pressure (atm) 

Ra agglomerate radius (m) 

Ronm specific ionic resistance (Q m) 

sS Laplace-transform variable (s = iw) 

Teel cell temperature (K) 

XO, molar fraction of oxygen 

Z impedance (Q m?) 

Greek letters 

€B backing porosity 

T ratio of the active surface area inside the agglom- 
erate to unit agglomerate volume (m7!) 

n overpotential (V) 

A product of the agglomerate surface area (4 R2) 
and the density of the active agglomerates in the 
catalyst layer (m7!) 

w angular frequency (s7!) 


assumed that the catalyst material in the fuel cell electrode is 
contained within cylindrical agglomerates [3]. The cathode cat- 
alyst layer contains agglomerated zones, filled with electrolyte. 
The agglomerates are intercalated with open hydrophobic chan- 
nels through which the oxidant reaches the surface of the flooded 
pores, dissolves, and diffuses to the catalyst particles, where the 
Faraday step occurs. Springer and Raistrick expanded this model 
for the simulation of ac impedance spectra [4]. However, SEM 
plots indicate that the assumption of spherically shaped agglom- 
erates is more realistic for PEM fuel cell electrodes [5]. For 
simulating polarisation curves, agglomerates of spherical shape 


in the cathode electrode are assumed to be present in [6-8]. 
Jaouen and Lindbergh upgraded their steady-state agglomerate 
model to simulate ac impedance spectra [9]. For simplification, 
they neglected the gas diffusion limitation in the backing layer. 

The mechanism of the oxygen reduction reaction in the fuel 
cell cathode is discussed by a number of authors. A doubling 
or even fourfold increase of the apparent Tafel slope has been 
reported and controversial interpreted. One explanation could 
be a change of the reaction mechanism. Damjanovic proposed 
the formation of several intermediate species during the oxygen 
reduction reaction (ORR) [10-12]. It is also well known that 
mass-transport limitations in the fuel cell cathode can cause a 
change in the apparent Tafel slope. When the cathode is con- 
trolled by the ORR and the diffusion of oxygen in the catalyst 
layer or by the ORR and the proton migration in the catalyst 
layer a doubling of the Tafel slope is observed [7]. An alternative 
explanation can be given by assuming spherical agglomerates in 
the electrode. According to this approach, the change of the Tafel 
slope is not caused by diffusion limitation of oxygen through 
the catalyst layer but by diffusion limitation of oxygen into the 
agglomerates [5,7,8]. 

In this paper, we present model-based analysis to extract 
physical and electrochemical parameters of the fuel cell cathode 
from impedance measurements. We developed a new cathode 
model based on the assumption of spherically shaped agglom- 
erates. The model based on a coupled differential equation 
system describing: ternary gas diffusion in the backing (O2, No, 
water vapour), Fickian diffusion and Tafel kinetics for the oxy- 
gen reduction reaction (ORR) inside the agglomerates, proton 
migration with ohmic losses and double-layer charging through 
the catalyst layer. The interaction between the electrode and 
the backing layer is taken into account for the calculation of 
ac impedance spectra and steady-state polarisation curves. Our 
model represents mass-transport limitation in the agglomerates 
as well as in the backing layer. The latter is significant at medium 
and high current density values and essential for simulating the 
behaviour of a PEM fuel cell in the whole current density range. 
According to our model, the change in the apparent Tafel slope 
is due to diffusion limitation of oxygen inside the agglomerates. 

The model is used for a detailed analysis of mass-transport 
losses and losses due to the slow ORR in the cathode of a com- 
mercially available membrane electrode assembly from W.L. 
Gore & Associates (PRIMEA® Series 5510). We extract physi- 
cal and electrochemical parameters by comparing measurement 
data and modelling results. Moreover, our cathode model can be 
used in more comprehensive mathematical models that include 
the fluid dynamics in the flow-fields. 


2. Modelling 


Our model takes into account both the structure of the elec- 
trode and the most important physical and electrochemical 
phenomena that occur in the electrode. These are the reaction 
kinetics of the ORR and the mass-transport limitation of oxygen 
caused by diffusion in the agglomerates. Furthermore, the model 
accounts for multicomponent gas diffusion in the backing layer. 
In particular, the interaction of the gas diffusion layer (GDL) 
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Catalyst Layer PEM 


Agglomerate 


Fig. 1. Schematic diagram of the spherical agglomerate cathode model. 


and the porous gas diffusion electrode of state-of-the-art PEM 
fuel cells is addressed. 

The model considers only the cathode side of a PEMFC, 
losses on the anode side are neglected because of the small polar- 
isation of a hydrogen electrode. The inner structure of the porous 
gas diffusion electrode is taken into account. In the catalyst layer, 
the dispersed catalyst (Pt/C) forms agglomerated zones filled 
with the electrolyte. These zones are surrounded by hydrophobic 
channels through which the reactant gases can diffuse. Agglom- 
erates of mean radius Ra are assumed to be homogeneously 
distributed within the catalyst layer. Fig. 1 shows a schematic 
diagram of the cathode. A position-dependent activation overpo- 
tential is accounted for. The centre of an agglomerate determines 
its position within the catalyst layer. The size of the agglomer- 
ates is sufficiently small compared to the thickness of the active 
layer, so that an uniform overpotential within each agglomerate 
can be assumed. 

The reactant gas dissolves in the bounded water of the 
ionomer at the boundary between the open gas pores and the 
agglomerates. Then the oxygen diffuses through the humidified 
electrolyte inside the agglomerates to the active sites where the 
electrochemical reaction takes place. The diffusion coefficient 
of the dissolved gas in the humidified electrolyte is about three 
orders of magnitude smaller than the diffusion coefficient of 
gas in the open pores. Hence, we assume that diffusion limita- 
tion in the open gas pores is negligible compared to diffusion 
losses within the agglomerates. Therefore, the oxygen concen- 
tration outside the agglomerates is set constant over the whole 
catalyst layer. Oxygen depletion occurs only inside the agglom- 
erates where the electrochemical reaction takes place and in the 
backing layer due to mass-transport limitation. 

At the moment this model does not describes the effects of 
liquid water. Thus, we are not able to predict flooding of the 
porous media or dehumidification of the ionomer. But these 
mentioned loss mechanisms can be interpreted by a change of 
model parameters like the effective porosity or the ionic resis- 
tance. 


2.1. Gas transport in the backing layer 


The model of the backing layer is based on the work of 
Springer et al. [2] where the equations are discussed in more 
detail. For better understanding, we review the most important 


assumptions and equations governing the gas transport in the 
backing layer. 

Ternary diffusion of oxygen, nitrogen and water vapour 
occurs through the backing layer. Interdiffusion of these gases 
through the porous backing is calculated by means of the 
Stefan—Maxwell equation. Eq. (1) relates the instantaneous gra- 
dient of the oxygen mole fraction to the fluxes of oxygen, 
nitrogen and water vapour, 


dxo, 
dy = [æo (x0, Nw aR XwNo,) JF Xo, NNa 
a )No,] : (1) 
= (1 = xo, — Xw) ol] -Rn > 
: i cDeo,N, 
where c and ao, are given by 
1 Py T D 
=, ato, = R, (2) 
Vm Ps Teell Dsoow 


A Bruggeman expression for Deo,n, includes the effective 
porosity €g and the tortuosity t of the gas diffusion layer: 


1.823 
Teen 
Ts i 


DeoN, = 2 bso ( (3) 
Here, x; is the mole fraction of the gas species i, N; is the molar 
flux of the gas i, Dsij is the binary diffusion coefficient at stan- 
dard conditions, 71 is the cell temperature, Ts = 273.16 K is 
the standard temperature, vm = 22.4 x 1073 m? mol! is the 
molar gas volume, c is the total gas concentration in the backing, 
P, is the inlet pressure and Ps = 1 atm is the standard pressure. 
Equilibrium saturation of the cathode air stream with water 
vapour and isothermal conditions are assumed within the back- 
ing layer. Moreover it is assumed that the total pressure, i.e. the 
sum of the partial pressure of oxygen, nitrogen and water vapour 
remains constant inside the backing layer. Due to this constant 
total pressure within the backing layer the mole fraction of the 
water vapour is constant. Therefore, the gradient of water vapour 
in the backing layer vanishes in the Stefan—Maxwell equation: 


dxw 
‘dy = [xo (xwNo, = XO, Nw) F QN, (Xw NN, 
ad )N ie 0 (4) 
— (1 — xo, — x =0, 
Oo w)tVw Dios 

where 

DsoyN> 
en, = =, 

Dswn> 


With an instantaneous ionic current density jeen, there is an uni- 
form molecular flux (jce/4F) of oxygen and nitrogen through 
the backing layer. Thus the nitrogen flux can be replaced by the 
difference between jce /4F and the oxygen flux. With this sub- 
stitution in Eq. (4), we can determine the water vapour flux Ny 
which is dragged along with the oxygen flux: 


Xw Ny, + xw QO, /AN, Cee /4F — Ny) 


Ny = 
1 — xy + (@0,/an, — 1) xo, 


(5) 
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Inserting this water vapour flux in Eq. (1), we obtain an expres- 
sion for the nitrogen flux: 


cDeo.n, | | jeen dxo, 


Nn, = lgG Z|, (6 
N2 b Ala (B1 — x0) + lB G(xo,) T (6) 
where 
Bi + 2x0, 
G ae. = xy, 
(xo) EEG By Xy 
a 
po = > —1, B3=1—xw + 00,%w, 
QN, 


and the characteristic backing current density is given by 
ig = 4 FcDeo,n, / lB. The nitrogen flux throughout the back- 
ing layer Ny, given by Eq. (6) is assumed to vanish under 
steady-state conditions because there is no nitrogen flux at the 
backing/catalyst layer interface. This leads to an implicit equa- 
tion for the oxygen mole fraction at the interface backing/catalyst 
layer xọo,(lg) as a function of the instantaneous ionic current 
density jce1 through the catalyst layer of the fuel cell: 


jeet By — 226) (£ + po zam) 
= l —— l — A ], 
poea Ẹ So) se, 0 


where 


_ Bid. + Bo) B3 — Bi 


MS ie b > Bi Bi + Bs 


Although in our experiments the relative humidity of air at the 
inlet was clear under 100% (see below) the assumption of sat- 
urated air was made due to receive this implicit equation (Eq. 
(7)) for the oxygen mole fraction at the interface backing/catalyst 
layer and therefore a fast solution. Numerical calculations of the 
Stefan—Maxwell equation without these simplifications show 
that the error in the oxygen mole fraction at the interface is 
small by using the analytical solution given in Eq. (7). More- 
over, the water generation of the ORR humidify the gas in the 
backing. Therefore, the relative humidity of air is not constant 
but a function of the current density. This feature and the related 
characteristics like water condensation or liquid water transport 
does not appear in the model. 


2.2. The catalyst layer agglomerate model 


The oxygen flux inside an agglomerate is modelled using 
Fick’s first law: 


x9 CY, y, t) 
or æ ’ 


jb, (r, y, ) = —De (8) 
where jo, is the oxygen flux inside a flooded agglomerate, D is 
the diffusion constant of dissolved oxygen within an agglomer- 
ate and c/c* is the normalised local oxygen concentration with 
reference to the dissolved oxygen concentration on the agglom- 
erate surface c*. The oxygen mole fraction dissolved in water 
bounded in the ionomer, and consequently the dissolved oxygen 
concentration c*, is calculated from the oxygen partial pressure 


Po, in the open gas pores using Henry’s law: 

A Pix 

diss _ PO2 a+Q, 
= — = 5 9 

XO, K K (9) 
where K = 4.34 x 104 atm is the Henry constant [13]. The oxy- 
gen partial pressure can be calculated with the oxygen mole 
fraction at the backing/catalyst interface and the inlet pressure 
from the backing equations of Springer et al. [2]. 

The mass balance equation combines the local change of 
oxygen flux with oxygen consumption due to the ORR and the 
change of oxygen concentration with time: 


Lə ag dg Po ( ORD 2.3 n(y, t) 
a ap Joy, N= nF of exp a 
d (clr, y,t) 
-o3 ( 2 ) l (10) 


where r is the radial coordinate within an agglomerate. The vol- 
umetric exchange current density iY = io I is the product of 
the exchange current density ig of the ORR at an Op partial 
pressure of | atm and the ratio I” of the active surface area 
inside an agglomerate to unit agglomerate volume. It is dif- 
ficult to extract J” from measurements, so we use in as a fit 
parameter. po, is the oxygen partial pressure in the open gas 
pores of the catalyst layer, 7 is the overpotential and b is the 
Tafel slope. The boundary conditions are (c(Ra, y, f)/c*) = 1 
and d/dr((c(0, y, t))/c*) = 0. Eqs. (8) and (10) determine the 
distribution of the normalised oxygen concentration inside the 
agglomerates depending on the local overpotential in the catalyst 
layer. 

The proton flux inside an agglomerate Jp is given by the 
charge balance equation: 


l ð iT . c(r, y, t) 2.3 ny, t) 
Dal Ilr y, D= -po ( a ) exp ( . 


b 
(1) 


The boundary condition for Eq. (11) is that the proton flux in 
the centre of an agglomerate is zero: A y, t) = 0. 

At steady-state, the differential equations (8)—(11) yield the 
current density on the surface of an agglomerate Jp(Ra, y)asa 
function of the position in the catalyst layer: 


f nFDc* f 2.3 n(y) 
jp(Ra, y= R yn PO, Dc* exp ( 5 
a 


VEX po. exp(2.3n(y)/b)Ra 
Vn FDc* 


The charge balance equation along the y-direction is given by 


x coth (12) 


djp(y, t) INY, t) 
dy ot 


where the ionic current density at the surface area of an agglom- 
erate is given by Jp(Ra, y, t), and the agglomerate centre is 
located at position y in the electrode. CpL,y is the volumetric 
double-layer capacitance of the electrode and A is the product 


= jp(Ra, y, NA + CoL,v (13) 
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of the agglomerate surface (4 R2) and the density of the active 
agglomerates in the catalyst layer. 

Ohm’s law is used to calculate the overpotential distribution 
in the catalyst layer: 


any, t) 
oy 


= Romm 2.0), (14) 


where Ronn is the specific ionic resistance of the catalyst layer. 
Due to the high electric conductivity of carbon compared to the 
low ionic conductivity of the ionomer, the ohmic losses in the 
carbon support are neglected. The proton flux at the interface 
between the backing layer and the active layer vanishes, leading 
to the boundary condition B (0, t) = 0. The second boundary 
condition ip (L, t) = jeen implies the initial condition for n(0, t) 
for the integration of the differential equation system. The solu- 
tion method is described in the next section. 


2.3. Laplace-transformation of the governing equations 


To calculate an impedance spectrum, each solution variable 
of the system of the differential equations (Eqs. (8)—-(14)) is 
perturbed about a steady-state solution. For a small perturba- 
tion, the response of the system can be deduced from the Taylor 
expressions limited to the first order. Introducing these Taylor 
expressions in Eqs. (8)—(14) and subtracting the steady-state 
equations, we obtain a Laplace-transformed system of differ- 
ential equations depending on frequency. The time-dependent 
perturbed variables are denoted with a bar, variables without 
a bar are the steady-state solutions and s = iw is the complex 
frequency. 

According to Springer et al. [2], an approximate analytical 
transfer function (Eq. (15)) for the dynamic response of the oxy- 
gen mole fraction xo, (y = lg, s) at the interface between the 
backing and the catalyst layer can be found depending on the 
dynamic response of the ionic current density jce(s) and the 
steady-state oxygen mole fraction te: This equation is used 
for the coupling of the backing with the catalyst layer in the 
frequency domain 


Jeeil(S) Xo, (lp) — Bi 
ig  G(xo,(lB)) 


V5 b1 B/ G0, (B) eons 
V5 b1 / G0; (lB) Deon 


X07 (lB, s) = 


x | tanh 


(15) 


Laplace-transformation of Fick’s first law describing the dif- 
fusion of dissolved oxygen into an agglomerate yields: 


cr, y, 8)/c*)  — elr, y)/ 2) 
c +c 5 
or or 


ace o( 


(16) 


where c(r, y, s)/c* denotes the dynamic response of the nor- 
malised local oxygen concentration. 


The Laplace-transformed mass balance equation reads: 


1 30? 76, y, 5) 
r2 or 


ae cn S)\ oyp (2310) 

~ nF” c* p b 
iY — fcr y) a 2.3 n(y) 
nF” c* j b 
iY clr, y) \ 2.3 2.3 n(y) 
nF’ c* Be b 


sunen h, (17) 


c* 


where the dynamic response of the oxygen partial pressure is 
given by Po, = Xo, (lB, $) Pa. 

For the charge balance equation inside an agglomerate we 
obtain: 


1 Ar? 7507, y, 8) 
r2 or 


— eV c(r, y, s) 2.3 n(y) 
= —ip Po2| — — } xP | -z 


v___¢(r, y) oth (= o) 


to POs c* b 


y c(r, y) 2.3 2.3 n(y) 
iù PO» a r exp 5 


) NY, $). (18) 
The boundary conditions for Eqs. (16)—(18) are as follows: 


e there is no variation of the dynamic response of the normalised 
oxygen concentration on the surface of the agglomerate 
(c(Ra, y, 8)/c*) = 0, 

e the gradient of the dynamic response of the normalised oxy- 
gen concentration in the centre of the agglomerates vanishes 
a(c(O, y, s)/c*)/or = 0, 

e the dynamic response of the current density vanishes in the 
centre of an agglomerate, i.e. EO. y, t)=0. 


Inserting Eq. (16) into Eq. (17) and solving this differen- 


tial equation, we obtain (2) as a function of 7(y, s) and 


n(y). Inserting this expression for (2) into Eq. (18), we 
obtain an analytical solution for the dynamic response of the 
current generation per agglomerate F5(Ra, y, s). The dynamic 
current response of the agglomerates is used to calculate the 
impedance of the entire cathode layer. Therefore, the differen- 
tial equation system containing the Laplace-transformed charge 
balance equation in the y-direction (Eq. (19)) and Ohm’s law 
(Eq. (20)) has to be solved. 


aO) _ 


ay F5(Ra, y, s)A F SsCoL vad, s), (19) 
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an(y, s) 
dy 


= Rohm jp; 5). (20) 
The boundary conditions are 


° j (0, s) = 0: no variation of the current density at the back- 
ing/catalyst interface, 

e 7(0, s) = 1: due to the linearity of the differential equation 
system, an arbitrary number not equal to zero can be chosen 
for 7(0, s). 


The cathode impedance Z,atn(@) at a frequency of w is calcu- 
lated by the quotient 7(L, iw)/ ib (L, iw). For comparison with 
measured spectra the measured high frequency resistance Ryp 
is added to the cathode impedance 


Z(w) = Zeath(w) + Rpr. (21) 


We solve the differential equation system using 
MATHEMATICA 4®. The steady-state equations are solved 
as follows: the oxygen mole fraction at the backing/catalyst 
interface is determined by means of the implicit equation (Eq. 
(15)) for a given current density jcey. With the knowledge of 
the dissolved oxygen concentration on the agglomerate surface 
(Eq. (9)), we obtain an analytical expression of the current 
generation per agglomerate as a function of the overpotential 
n(y) (Eq. (12)). This term is used as a source term for the 
current density in the catalyst layer in Eq. (13). The differential 
equation system (Eqs. (13) and (14)) can be solved numerically 
using the function NDSolve. In doing so, the initial condition 
7(0) is adapted using the Newton method (FindRoot) in order 
to obtain the second boundary condition i (L) = jeen. The cell 
potential is calculated by 


Usen = 1.23V — nL). (22) 


NDSolve is a built-in function of MATHEMATICA® which 
finds numerical solutions of differential equations and sys- 
tems using an Adams Predictor—Corrector method for non-stiff 
differential equations and backward difference formulae (gear 
method) for stiff differential equations. It switches between the 
two methods using heuristics based on the adaptively selected 
step size [14,15]. The results are given in terms of interpolating 
functions. Thus, we obtain the cell potential for a given cur- 
rent density and the distribution of the overpotential within the 
catalyst layer. 

The catalyst layer equations in conjunction with the back- 
ing layer transfer function (Eq. (15)) are solved to calculate 
the cathode impedance Z.arn. NDSolve is used to solve the 
Laplace-transformed equations (Eqs. (19) and (20)). Because 
of the interaction between the backing layer and the catalyst 
layer, the value i (L, s) has to be adapted iteratively using the 
Newton method (FindRoot) so that a self-consistent solution to 
the backing transfer function is found, i.e. i (L, s) = Jeenls). 


3. Experimental 


A small test fuel cell was constructed to validate our agglom- 
erate model. To guarantee homogeneity in current density, 
temperature and gas supply over the whole MEA, the geomet- 
ric active area is only 1.4cm x1.4cm. The cell temperature 
can be controlled by water channels in the graphite plates that 
contain the flow-field. Serpentine flow-fields were machined in 
the graphite plates on both the anode and the cathode side. 
An untreated gas diffusion layers of type TORAY®TGP — 
H — 120 were utilised. A commercial MEA from GORE™ 
(PRIMEA® Series 5510, membrane thickness 25 um) was 
investigated. The cell was connected to a Solartron 1286 elec- 
trochemical interface in combination with a Solartron 1255 HF 
frequency response analyser. The impedance spectra were mea- 
sured in the frequency range between 0.1 Hz and 10 kHz. All ac 
impedance spectra were measured between the cathode and the 
anode, i.e. without a reference electrode. The gas flow was regu- 
lated with a gas flow controller, with the flow rate set to 110 sccm 
for hydrogen and 200 sccm for air at atmospheric pressure. Dry 
hydrogen was used on the anode side, on the cathode side the 
air was humidified with a wash bottle to 100% relative humid- 
ity (RH) at room temperature (20°C). The measurements were 
made at a cell temperature of 39°C. At this temperature the 
cathode gas has a RH of 33.4%. 

The cell potentials of all polarisation curves mentioned in 
the text are corrected for the ohmic drop in the membrane. To 
determine the resistance of the membrane we used the real part of 
the 10 kHz impedance. For the parameter extraction all observed 
currents were corrected for an estimated crossover current of 
7Am?[16]. 


4. Comparison of measurements and simulation 


The dotted lines in Fig. 3 shows measured impedance spectra 
in the potential range between 0.75 and 0.5 V. The high- 
frequency resistance decreases with decreasing cell potential. 
This can be explained by a change of the MEA water content. 
At low cell potential, i.e. with higher current densities, the prod- 
uct water humidifies the MEA leading to a smaller protonic 
resistance. 

To validate our model, we searched for a set of parameters 
such that the model reproduces both measured impedance spec- 
tra at different current density values and the experimentally 
observed polarisation curve. The following parameters were 
fixed during the fitting procedure: 


e The thickness of the catalyst layer was extracted from a micro- 
scope image to be le = 10 pm. 

e The diffusion coefficient of dissolved oxygen in an agglom- 
erate D is set to 2 x 107° m? s71, according to Springer et al. 
[4] 

e A gas diffusion layer with a thickness of 370 pm was used 
in the experiment. It is assumed that the thickness of the gas 
diffusion layer is reduced to [g = 350 wm due to the contact 
pressure applied to the cell, thus the latter value was used for 
the simulation. 
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Fig. 2. Measured and simulated polarisation curves for low current density. 
The simulation shows good agreement with the measurement. A change of the 
apparent Tafel slope is obtained. 


e The value of the agglomerate radius (Rg = 0.1 wm) was 
adopted from Ihonen et al. [5]. 

e A value of 6.5 is chosen for the backing tortuosity t (Springer 
et al. [2]t = 7). 

e During the measurement, the operating pressure P, was 1 atm 
(pressure drop in the flow-field is negligible) and the cell 
temperature was controlled to be 312 K. 


Beginning with a comparison of the measured and simu- 
lated polarisation curves for low current density values results 
in a Tafel slope b of 78.3mV and a volumetric exchange 
current density iy of 7.74 x 10° Am™? (see Fig. 2). Mass- 
transport limitations of dissolved oxygen into the agglomerates 
start already at low current density values indicated by the 
observed doubling of the apparent Tafel slope. The polarisa- 
tion curve show also a sudden drop of the cell potential in 
the high current density range. This mass-transport limitation 
starting at the cell potential of 550 mV is attributed to flood- 
ing effects in the GDL and CL. For the following calculations 
these two parameter values were hold constant. We fitted the 
computed impedance spectra to the measured spectra in the 
range of low current density values (639 Am?, 1388 Am’), 
where the spectra are not affected by oxygen diffusion limita- 
tion in the backing layer. Therefore, the backing layer parameters 
have no influence on the spectra under these operation condi- 
tions. 

The fit parameters are A, Rohm and CpL,y. Each parameter 
reveals a characteristic influence on the frequency-dependent 
shape of the simulated impedance spectrum. Moreover, the 
parameter influence on a simulated impedance spectrum changes 
depending on the overpotential. Therefore, the values of the 
parameters can be extracted by simultaneously fitting several 
simulated spectra to measured impedance spectra that corre- 
spond to different overpotential values. We obtain the best 
fit with a value of A of 35,000m7—!. For further fits with 
higher current density values, where the spectra are affected 
by the diffusion limitation of the backing, we fixed the value 
of A and added the backing porosity eg as fit parameter to 
Ronm, CDL,v- 
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Fig. 3. Comparison between simulated (lines) and measured (dots) spectra for 
different current densities from 639 A m~? to 4612 Am7?. 


Fig. 3 shows the simulated and measured spectra for dif- 
ferent current density values from 639 A m7? to 4612 Am7?. 
Experimentally, the diameter of the impedance arc decreases 
with increasing current density due to a reduced charge transfer 
resistance for higher overvoltages. At a current density value of 
2600 A m7? the measurement reveals a minimum value of the 
radius of the impedance arc, indicating that the charge trans- 
fer resistance is not longer the dominating loss mechanism. 
Subsequently the radius grows again as the current density 
increases. Oxygen mass-transport limitations are the main loss 
in this current density range. This characteristic behaviour of 
the impedance arcs is reproduced by our model. The simulated 
impedance spectra show good agreement with the measured 
spectra in the high and medium frequency range. The simulated 
impedance values for low frequency are too small. The measured 
impedance spectra reveal an additional mass-transport resistance 
that is not accounted for by our model. For instance, this could 
be diffusion limitation of gaseous oxygen in the open gas pores 
or the anodic impedance [17]. 

Fig. 4 shows the normalised values of the measured high fre- 
quency resistance Ryp and the normalised fit results of Rohm 
and CpL,y for different current density values. The maximum 
values used for the normalisation can be found in Table 1. It can 
be clearly seen that the measured high frequency resistance Rup, 
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Fig. 4. Normalised values of the measured high frequency resistance Ryp and 
the normalised fit results of Ronm and CpL,y for different current density values. 
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Table 1 

Maximum value of the high frequency resistance (HFR), the catalyst layer resis- 
tance (Ronm) and the double layer capacity (Cpi,v) used for the normalisation 
in Fig. 4. The HFR is experimentally measured and Ronm as well as CpL,v are 
obtained from the fit procedure. 


High frequency resistance (x 10~* Q m?) 5.6 
Catalyst layer resistance (Q m) 5.49 
Double layer capacity (x 107 F m7?) 1.9 
Table 2 

Parameters obtained from the fit procedure 

b (V) 78.3 

iy (Am™°) 7.74 x 105 
A w7) 35,000 
Ronm (Qm) 2.02-5.49 


CpL.v (Fm™’) 6.8 x 106—1.9 x 107 


contributing mainly from the membrane, corresponds to the fit- 
ted ohmic resistance of the catalyst layer Ronm. For low current 
density, the water generation due to the ORR is not sufficient for 
a high stoichiometric air flow. This leads to a dehumidification 
of the MEA and thus to a decrease of the ionic conductivity in 
the catalyst layer and in the membrane. The humidification level 
of the MEA increases with increasing current density, leading to 
an decrease of Ronm and Rpr. Additionally, a better humidifi- 
cation of the electrode leads to a higher volumetric double-layer 
capacity CpL,v, indicating an increase of the three-phase bound- 
ary. 

At high current density the impedance spectra show no 
change in the high frequency resistance, indicating that the 
membrane is well humidified. Flooding occurs in the region of 
the limiting current density. The open pores become filled with 
liquid water so that the open pore space available for gas flow 
decreases. By setting the backing porosity €g to an value of 
13%, we are able to fit the spectra at high current density values. 
This unrealistic value of e shows that our one-phase model is not 
valid for values near the limiting current density if liquid water 
is present in the porous media. The parameter obtained from the 
curve fitting are summarized in Table 2. A fuel cell model that 
covers the whole range of current density values from zero to 
limiting current density without changing model parameters has 
to include the water management. The simulation of flooding 
effects and dehydration of the ionomer are essential for the 
interpretation of experimental data. The implementation of a 
transport equation for liquid water in the porous media (GDL, 
CL) coupled with water film formation around the agglomerates 
lowering the oxygen diffusion to the active sides and the descrip- 
tion of the water content in the ionomer which contribute to a 
change in the proton conductivity and activity [18] is ongoing 
work. 


5. Parameter study 


This section presents an investigation of the effect of the cath- 
ode parameters on the shape of the polarisation curves and the 
impedance spectra. The parameter study outlined here helps to 
interpret measurement results and to analyse the influence of the 


cathode parameters on the shape of the curves within different 
measurement regimes. 


5.1. Polarisation curves 


The parameter study was performed by varying each model 
parameter by about +50% with respect to the following baseline 
case: b = 53 mV, Ra = 0.1 ym, D=5 x 10719 m? s—!, iy = 
15.2 x 104 Am=3 atm7!, le = 10 um, A = 5 x 10° m7!. Due 
to the strong influence of the Tafel slope on the polarisation 
curve, this parameter was only varied by about +10% of the 
baseline case value. 

First the effect of the Tafel slope is studied. Fig. 5(a) shows 
simulated polarisation curves for three different values of the 
Tafel slope. As mentioned before, the curves show a strong 
response to the change of the Tafel slope over the whole cur- 
rent density range. A doubling of the apparent Tafel slope is 
observed in all cases. 

At low current density values, the performance of the cath- 
ode depends strongly on the agglomerate radius (Fig. 5b). Here, 
the ORR occurs within the entire volume of the agglomer- 
ate. Diffusion limitation does not occur here. An increase of 
the agglomerate radius causes an increase of the reaction vol- 
ume. The diffusion limitation is more pronounced for a higher 
agglomerate radius. The doubling of the Tafel slope can be 
observed at a lower current density in this case. 

Fig. 5(c) shows the dependence of the polarisation curve on 
the diffusion coefficient of dissolved oxygen into the agglomer- 
ates. For low current density, the oxygen diffusion inside the 
agglomerates is not a limiting effect. Thus, the polarisation 
curves corresponding to low current density do not depend on 
the diffusion coefficient. In this case, the whole agglomerate 
is supplied by oxygen and utilised for the ORR. The diffusion 
coefficient determines the current density value where the dif- 
fusion limitation occurs and the apparent Tafel slope changes. 
A smaller diffusion coefficient causes a shift to smaller values 
of the current range corresponding to the doubling of the Tafel 
slope. 

The exchange current density is varied in Fig. 5(d). An 
increase of the exchange current density shifts the polarisation 
curves to lower overpotential. No influence on the doubling of 
the apparent Tafel slope is observed. 

The product A of the agglomerate surface area and the density 
of the agglomerates and the thickness of the catalyst layer L have 
almost the same influence on the polarisation curves (Fig. 6). 
Both a thickness increase of the catalyst layer (Fig. 6a) and an 
increase of A (Fig. 6b) lead to a higher volume that is available 
for the ORR. Consequently, more current can be generated in 
the catalyst layer at the same overpotential. 

An influence of the backing layer parameters on the polar- 
isation curve is only evident at high current density. Here, 
oxygen diffusion through the backing is limiting the fuel cell 
performance. Accumulation of condensed liquid water deter- 
mines the saturation in the gas diffusion layer. As explained 
in Section 4, this effect can be interpreted as a decrease 
of the effective porosity of the gas diffusion layer in our 
model. 
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Fig. 5. Parameter variations show the influence on the doubling of the apparent Tafel slope. (a) Variation of the Tafel slope, (b) variation of the agglomerate radius, 
(c) variation of the diffusion coefficient, (d) variation of the exchange current density. 
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Fig. 6. Variation of the thickness of the catalyst layer (a) and of the ratio of the agglomerate surface area to unit volume (b) shows a similar influence. 


5.2. Impedance spectra 


Fig. 7 shows the effect of varying the ionic resistance in the 
catalyst layer. The spectra are calculated at a current density 
of 639 Am~?. The higher the ionic resistance of the catalyst 
layer, the more distinct is the linear branch with an angle of 
45° between the real axis and the impedance plot in the Nyquist 
representation. This branch reflects coupling of the distributed 
ionic resistance and the distributed capacitance of the catalyst 


layer. That is, the double-layer charging effects and the proton 
transport dominate the overall electrode response in the high- 
frequency regime [2,19]. 

Fig. 8(a) shows simulated impedance spectra for different 
current density values. The parameters of the baseline case are 
close to the validated parameter values (Table 2). For low cur- 
rent density values the spectrum is determined by the ORR, 
indicated by the decrease of the Nyquist arc radius with increas- 
ing overpotential. This behaviour changes for high current 
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Fig. 7. The influence of the ionic resistance in the catalyst layer on the impedance is shown. 
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Fig. 8. Investigation of the processes that dominates the shape of the impedance spectra. (a) Baseline case, (b) reduced: CpL,y, (c) reduced: CpL,y and ep, (d) 


reduced: D. 


density values where the radius of the Nyquist arc increases 
with increasing overpotential, indicating that the charge trans- 
fer resistance is not longer the dominating loss mechanism. 
To understand this change we reduced the volumetric double- 
layer capacity CpL,v from 107 Fm~? to 10° Fm~? in Fig. 8(b) 
without changing other parameters. Two separate arcs appear 
at higher current density values. The high-frequency arc corre- 
sponding to the ORR decreases with increasing overpotential. 
The low-frequency arc increases with increasing overpotential. 
This arc corresponds to the backing impedance. Fig. 8(c) shows 
the same spectra like in Fig. 8(b) but simulated with an addi- 
tional lower effective backing porosity, changed from 0.14 to 
0.12. The lowered backing porosity does not influence the high- 
frequency arc but the radius of the low-frequency arc is enlarged 
at high current density compared to the low-frequency arc in 
Fig. 8(b). In addition, Fig. 8(d) shows impedance spectra with a 
minimised diffusion coefficient for oxygen inside the agglomer- 
ates from 5 x 107° m’s~! to 5 x 10710 m? s7! (CpL,y has still 
the reduced value of 106 Fm~+). Now, the high-frequency arc 
expands due to diffusion limitation in the agglomerates. From 
the results of Fig. 8 it can be concluded that a high volumetric 
double-layer capacity masks the two distinct arcs (like in the case 
of our investigated MEA) and therefore the loss mechanism can- 
not separated easily. Depending on the current density either the 
first high frequency arc or the second low-frequency arc, which 
are merged to one arc in case of a high double-layer capac- 
ity, dominates the overall impedance spectrum and therefore 
the radius decreases or increases with increasing overpotential. 
This explains the shape of the measured impedance spectra at 
different operation points. 


5.3. Sensitivity analysis 


By means of parameter variation, the influence of the model 
parameters on the impedance spectra at different current density 


values was investigated. The behaviour of the impedance spec- 
tra once with low current density (639 A m~?) and once with 
high current density (4341 A m~?) were tested. The validated 
parameter set was selected as a baseline case. Each parameter 
was varied about +10% in comparison to the baseline value. 
The parameters of the gas diffusion layer have no influence on 
the impedance spectra at low current density. This statement is 
certainly only valid, as long as the parameter values of the gas 
diffusion layer (GDL) do not give rise to a limiting current in 
the low current density region. 

In the low current density range, the electrode dominates 
the shape of the impedance spectra: a reduction of the prod- 
uct A of the agglomerate surface area and the density of the 
active agglomerates in the catalyst layer, a reduction of the dif- 
fusion coefficient of dissolved oxygen in the agglomerates D or 
an increase of the agglomerate radius Ra cause an increase of 
the diameter of the impedance arc. At low current density val- 
ues, the variable ¢ dominates the Eqs. (17) and (18), which are 
linked to the diffusion constant and the agglomerate radius by 
Eq. (16). 

In the high current density range, the parameters of the gas 
diffusion layer dominate the shape of the impedance spectra. 
The backing parameters determine the current density at which 
diffusion limitation in the GDL occurs. Thus an increase of 
the thickness of the gas diffusion layer lg or the tortuosity t 
or a reduction of the backing porosity €g leads to an increase 
of the diameter of the impedance arc at high current density 
values. 

This effect can be explained by means of Eqs. (17) and 
(18). The dynamic response of the oxygen partial pressure po, 
dominates the sink terms of the differential equations at high 
current density. po, is linked to the transfer function for the 
dynamic response of the oxygen mole fraction (Eq. (5) in [2]) 
by the relationship Po, = XO, Pa and consequently to the model 
parameters of the GDL. 
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The Tafel slope b has a strong influence on the shape of the 
impedance spectra over the whole current density range. An 
increase of the Tafel slope causes an increase of the diameter of 
the impedance arc. 

This can also be explained by means of Eqs. (17) and (18). On 
the right-hand side of Eqs. (17) and (18), there is an exponential 
function with the Tafel slope in each addend. Thus, the Tafel 
slope always has a strong influence on the shape the impedance 
spectra, regardless of which Laplace-transformed function is 
dominating the equation. 

It was found that the volumetric exchange current density is 
does not affect the spectrum. Due to the characteristic of the dif- 
ferential equation system the exchange current density iy only 
affects the absolute value of the overpotential 7 but not the shape 
of its distribution in the catalyst layer for a given current density 
value. Due to the fact that only the product of i exp(2.3n/b) 
appears in Eqs. (17) and (18), a change in iy does not affect the 
impedance because the product remains unchanged. A mathe- 
matical explanation of this fact can also be found in Jaouen et 
al., Appendix B [9]. The volumetric exchange current density 
only affects the current-voltage characteristics and can therefore 
only be extracted on polarisation curves. 

A larger oxygen gas pressure at the inlet improves the reaction 
due to a higher dissolved oxygen concentration in the agglom- 
erates. Therefore, the diameter of the impedance arc decreases. 

Variation of the volumetric double-layer capacity CpL,y over 
a small range does not change the shape of the impedance 
spectrum in the Nyquist representation. In the Bode represen- 
tation, one recognises a phase shift of the peak of the absolute 
impedance values towards lower frequencies by increasing the 
capacity. For a very low double-layer capacity at high current 
density two separate arcs occur in the Nyquist plot (see Fig. 8). 

Inthe high-frequency range, the ionic resistance Ronm and the 
thickness of the catalyst layer l affect the 45° branch, and there- 
fore also the regime of the impedance spectra at low frequencies 
(see Fig. 7). 

There are additional parameters in the model, that influence 
the shape of the impedance spectra in a similar way: 


e The product A of the agglomerate surface area and the density 
of the active agglomerates in the catalyst layer, the diffusion 
coefficient of the dissolved oxygen in an agglomerate D and 
the agglomerate radius R, have the same influence on the 
simulated impedance spectra: 


Ao Ds R. 


e Likewise, the influence of the tortuosity t, porosity €g and 
the thickness of the gas diffusion layer lg can hardly be dis- 
tinguished: 


TO ER Ol]. 


Therefore, ex situ measurements are required to determine 
parameters like porosity and the mean radius of the agglomer- 
ates, e.g. porosimetry measurement of the GDL and STEM plots 
of the cathode. 


6. Conclusion 


A one-dimensional cathode agglomerate model for analysing 
impedance spectra was developed and validated with experi- 
mental data. The model predicts the measured doubling of the 
apparent Tafel slope, which occurs due to diffusion limitation 
of dissolved oxygen into the agglomerates. The model shows a 
good agreement in the characteristic behaviour of the measured 
impedance spectra at different current density values. In the low 
frequency range the model predicts a too low impedance in com- 
parison with measured spectra. An explanation for this could be 
that the anodic impedance or the diffusion limitation of gaseous 
oxygen in the catalyst layer cannot be neglected. 

The fitting results show that a better humidification of the 
electrode leads to a higher volumetric double layer capacity. 
The fitted ohmic resistance of the catalyst layer shows the same 
behaviour on the current density as the measured high frequency 
resistance due to product water humidifying the ionomer. 

The influence of the model parameters on the polarisation 
curve and the impedance spectrum were investigated. 

The implementation of a transport equation for liquid water 
in the porous media coupled with water film formation around 
the agglomerates lowering the oxygen diffusion to the active 
sides is ongoing work. 
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